*****************************************
*** 2. LATE PSYCHOLOGICAL with CONTROLS
*** Prepare adjusted p-values
use "$gsdRawTemp\qvalue_ate_con_1", clear
foreach num of numlist 2/5{
	append using "$gsdRawTemp\qvalue_ate_con_`num'"
}
*append using "$gsdRawTemp\qvalue_atex_con_1"
*append using "$gsdRawTemp\qvalue_atex_con_2"
***Bonferroni q-values
qqvalue p , method(bonferroni) qvalue(qvalue)
gen bonferroni_q=round(qvalue, .001)
****Benjamini Hochbergen q-values
* Collect the total number of p-values tested
quietly sum p
local totalpvals = r(N)
* Sort the p-values in ascending order and generate a variable that codes each p-value's rank
quietly gen int original_sorting_order = _n
quietly sort p
quietly gen int rank = _n if p!=.
* Set the initial counter to 1 
local qval = 1
* Generate the variable that will contain the BH (1995) q-values
gen bh95_qval = 1 if p!=.
* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.
while `qval' > 0 {
	* Generate value qr/M
	quietly gen fdr_temp = `qval'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= qr/M
	quietly gen reject_temp = (fdr_temp>=p) if fdr_temp~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	quietly gen reject_rank = reject_temp*rank
	* Record the rank of the largest p-value that meets above condition
	quietly egen total_rejected = max(reject_rank)
	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
	replace bh95_qval = `qval' if rank <= total_rejected & rank~=.
	* Reduce q by 0.001 and repeat loop
	quietly drop fdr_temp reject_temp reject_rank total_rejected
	local qval = `qval' - .001
}	
quietly sort original_sorting_order

preserve
keep if family == "behavioral_controls_t1"
mkmat bh95_qval, matrix(B3) rownames(parm)
restore
keep if family == "behavioral_controls_t2"
mkmat bh95_qval, matrix(B4) rownames(parm)


************************************
*** 3. LATE PSYCHOLOGICAL with GEOGRAPHY CONTROLS
*** Prepare adjusted p-values
use "$gsdRawTemp\qvalue_ate_geocon_1", clear
foreach num of numlist 2/5{
	append using "$gsdRawTemp\qvalue_ate_geocon_`num'"
}
*append using "$gsdRawTemp\qvalue_atex_geocon_1"
*append using "$gsdRawTemp\qvalue_atex_geocon_2"
***Bonferroni q-values
qqvalue p , method(bonferroni) qvalue(qvalue)
gen bonferroni_q=round(qvalue, .001)
****Benjamini Hochbergen q-values
* Collect the total number of p-values tested
quietly sum p
local totalpvals = r(N)
* Sort the p-values in ascending order and generate a variable that codes each p-value's rank
quietly gen int original_sorting_order = _n
quietly sort p
quietly gen int rank = _n if p!=.
* Set the initial counter to 1 
local qval = 1
* Generate the variable that will contain the BH (1995) q-values
gen bh95_qval = 1 if p!=.
* Set up a loop that begins by checking which hypotheses are rejected at q = 1.000, then checks which hypotheses are rejected at q = 0.999, then checks which hypotheses are rejected at q = 0.998, etc.  The loop ends by checking which hypotheses are rejected at q = 0.001.
while `qval' > 0 {
	* Generate value qr/M
	quietly gen fdr_temp = `qval'*rank/`totalpvals'
	* Generate binary variable checking condition p(r) <= qr/M
	quietly gen reject_temp = (fdr_temp>=p) if fdr_temp~=.
	* Generate variable containing p-value ranks for all p-values that meet above condition
	quietly gen reject_rank = reject_temp*rank
	* Record the rank of the largest p-value that meets above condition
	quietly egen total_rejected = max(reject_rank)
	* A p-value has been rejected at level q if its rank is less than or equal to the rank of the max p-value that meets the above condition
	replace bh95_qval = `qval' if rank <= total_rejected & rank~=.
	* Reduce q by 0.001 and repeat loop
	quietly drop fdr_temp reject_temp reject_rank total_rejected
	local qval = `qval' - .001
}	
quietly sort original_sorting_order

preserve
keep if family == "behavioral_controls_t1"
mkmat bh95_qval, matrix(B5) rownames(parm)
restore
keep if family == "behavioral_controls_t2"
mkmat bh95_qval, matrix(B6) rownames(parm)
